Getting started

Adapted from: SOSTOOLS' SOSDEMO1 (See Section 4.1 of SOSTOOLS User's Manual) and Example 2.4 of [PJ08]

P. Parrilo and A. Jadbabaie Approximation of the joint spectral radius using sum of squares. Linear Algebra and its Applications, Elsevier (2008), 428, 2385-2402

using DynamicPolynomials
@polyvar x y
p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4

\[ 5y^{4} - x^{2}y^{2} + 2x^{3}y + 2x^{4} \]

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using SumOfSquares
import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver)
con_ref = @constraint(model, p >= 0)
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We see above that the solver found a feasible solution. We now inspect this solution:

q = gram_matrix(con_ref)

GramMatrix with row/column basis:
 SubBasis{Monomial}([y^2, x*y, x^2])
And entries in a 3×3 SymMatrix{Float64}:
  5.0                     4.4408920985006264e-17  -2.9518518518518517
  4.4408920985006264e-17  4.903703703703703        1.0
 -2.9518518518518517      1.0                      2.0

We can get the SOS decomposition from the gram matrix as follows:

sosdec = SOSDecomposition(q)

(-2.1236350925148186·y^2 + 0.685602733872361·x*y + 1.404328747693391·x^2)^2 + (-0.6931517747987985·y^2 - 2.1053486943394364·x*y - 0.020343251432301153·x^2)^2 + (0.0985627258797928·y^2 - 0.034050994900062616·x*y + 0.16567112157245356·x^2)^2

We now seek for the SOS decomposition of the following polynomial:

p = 4*x^4*y^6 + x^2 - x*y^2 + y^2

\[ y^{2} + x^{2} - xy^{2} + 4x^{4}y^{6} \]

We build the same model as previously with this new polynomial. Here we can use Model instead of SOSModel as we explicitly constrain p to belong to the SOS cone with p in SOSCone().

model = Model(solver)
con_ref = @constraint(model, p in SOSCone())
optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We can query the SOS decomposition directly from the constraint reference as follows:

sos_decomposition(con_ref)

(-0.8633544111907671·y + 0.19297758710767499·x*y + 3.0514841249382645e-18·x*y^2 + 1.9748040632461146·x^2*y^3)^2 + (-1.8744116508952374e-16·y + 0.7962855107812752·x - 1.0528501085696115e-16·x*y - 1.7452642425911875·x*y^2 - 8.635065971616514e-17·x^2*y^3)^2 + (0.24383463418983628·y - 1.7639164207923265e-16·x - 1.548639212897575·x*y + 3.741600489401518e-17·x*y^2 + 0.25793362243589013·x^2*y^3)^2 + (-5.3930412033044254e-17·y + 0.6049209744419549·x - 7.96731030959412e-17·x*y + 0.27599821010522496·x*y^2 - 9.028384658229797e-18·x^2*y^3)^2 + (-0.4417735074073057·y - 1.0850742441349411e-16·x - 0.10009637589813172·x*y + 1.468808219015887e-17·x*y^2 - 0.18335527863613862·x^2*y^3)^2

This page was generated using Literate.jl.